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5_^ Abstract. We extend the finite element method introduced by Lakkis and 

Oi , Pryer [2011] to approximate the solution of second order elliptic problems in 

^^ nonvariational form to incorporate the discontinuous Galerkin (DG) frame- 

^^ work. This is done by viewing the NVFEM as a mixed method whereby the 

QQ "finite element Hessian" is an auxiliary variable in the formulation. Represent- 
ing the finite element Hessian in a discontinuous setting yields a linear system 

I I of the same size and having the same sparsity pattern of the compact DG 

■^^ methods for variational elliptic problems. Furthermore, the system matrix is 

\—^ very easy to assemble, Thus this approach greatly reduces the computational 

^H complexity of the discretisation compared to the continuous approach. 

(^ We conduct a stability and consistency analysis making use of the unified 

■4-J framework set out in Arnold et. al. [2001]. We also give an apriori analysis 

y^ . of the method. The analysis applies to any consistent representation of the 

r^ finite element Hessian, thus is applicable to the previous works making use of 

I I continuous Galerkin approximations. 



\!7 1- Introduction 

*^i Nonvariational partial differential equations (PDEs) are those which are given 

j~^ in the form 

(N (1.1) -A:D2m = /, 

^T where X.'Y — trace(-X"'"V) is the Frobenious inner product between matrices. If 

the matrix A is differentiable then there is an equivalence between this problem 
and its variational sibling 

(1.2) - div(AVu) + DAVu = /, 



O 



where 



d d 



>< 
c5 (1.3) DA = f^a,a,,i(a;),...,^9,a,,d(a;)] . 

\i=l i=l / 

Rewriting in this form is sometimes undesirable. For example, if the coefficient ma- 
trix A has near singular derivatives the problem will become advection dominated 
and possibly unstable for conforming finite element methods. There is a wealth of 
material on the treatment of advection dominated problems |EG041 IESW051 c.f.]. 
If A is not differentiable then the problem has no variational structure. In this case 
standard finite element methods cannot be applied. 

In a previous work 'LPllj a finite element method for the approximation of the 
nonvariational problem ( |1.1[ ) was introduced. This involved the introduction of a 
finite element Hessian represented in the same finite element space as the solution 
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(modulo boundary conditions). The applications of the discrete representation of a 
Hessian of a piecewise function are becoming broader, for example, it can be used 
to drive anisotropic adaptive algorithms JAV021 IVMD+071 , as a notion of discrete 
convexity [AM09^ and in the design of finite element methods for nonlinear fourth 
order problems |Pryl2| . We are particularly interested in nonvariational PDEs due 
to their relation to general fully nonlinear PDEs 

(1.4) ^(D^u) = 0, 

which are of significant current research. There have been finite element meth- 
ods presented for this general class of problem for example in |Boh08) the author 
presents a C"'^ finite element method shows stability and consistency (hence con- 
vergence) of the scheme which requires a high degree of smoothness on the exact 
solution. In |FN09bl IFNOQaj the authors give a method in which they approximate 
the general second order fully nonlinear PDE by a sequence of fourth order quasi- 
linear PDEs. This is reminiscent of the vanishing viscosity method introduced for 
classically studying first order fully nonlinear PDEs. Efficiency of any method used 
to approximate a problem such as this is key. Each of the methods are computa- 
tionally costly due to their reliance on C^ finite elements |B6h08l [FNOQbj or mixed 
meth ods [F^ 9aj. 

In |BS91| a generic framework was set up to prove convergence of numerical ap- 
proximations to the solutions of fully nonlinear PDEs. This involved constructing 
monotone sequences of approximations which are typically applied to finite differ- 
ence approximations of the nonlinear problem [Obe06[ c.f.]. The assumption of 
consistency made in the |BS91j framework are incompatible with finite element 
methods, however, an extremely important observation made in |JS11| is that the 
consistency condition may be weakened to incorporate the finite element case using 
a localisation argument (in the case of isotropic diffusion) . 

In this contribution we present a method for the discontinuous approximation of 
the linear nonvariational problem ( |1.1[ ). We also present convergence analysis for 
a certain subclass of the nonvariational problems, those which are coercive. This 
allows us to use variational techniques to analyse the problem. We prove optimal 
convergence rates for the finite element solution in broken Sobolev norms. Note 
that the results presented here are immediately applicable to the method derived 
for the continuous case given in [LPllj . 

The algebraic formulation of the continuous approximation of the nonvariational 
problem requires the solution of large sparse (d+l) x N'^ linear system |LP1H 
Lem 3.3], where d is the dimension of the problem and N the number of degrees of 
freedom. Equivalently, using a Schiir complement argument, this can be reduced 
to an N'^ full linear system. The reason that this system is full is due to the 
global nature of the L2(ri) projection operator into a continuous finite element 
space. The motivation for extending the nonvariational finite element method into 
the discontinuous setting is the massive gain in computational efficiency over the 
continuous case. Indeed, due to the local representation of the projection operators 
in these discontinuous spaces we are able to make massive computational savings, 
in that the system matrix will become sparse and is the same size as that of a 
standard discontinuous Galerkin stiffness matrix. 

To test the method numerically we make use of the finite element package Dune 
[BBD+08al iBBD+OSb] . In this work we are interested in the asymptotic behaviour 
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of the discontinuous approximation. In a subsequent work we will study the com- 
putational gains using the discontinuous framework presented over the continuous 
one given in |LPllj . as well as exploit the powerful parallelisation capabilities of 
the package. 

The rest of the paper is set out as follows: In SJ2] we formally introduce the 
model problem and give a brief review of known classical facts about nonvariational 
PDEs. In ^ we examine the discretisation of the nonvariational method in the 
discontinuous Galerkin framework, making use of the unified framework set out 
in |ABCMd2] to derive a very general formulation of the finite element Hessian 
represented as a discontinuous object. We present some examples and examine 
the natural question of what happens when we try to eliminate the finite element 
Hessian from the formulation. In 34] we look at the consistency and stability of 
the finite element Hessian and present our main analytical results of convergence. 
Finally, in f|5] we detail a summary of extensive numerical experiments aimed at 
examining convergence and computational speed of the method presented. 

2. Problem formulation 

In this section we formulate the model problem, fix notation and give some 
basic assumptions. In addition we review the existence and uniqueness of the 
nonvariational problems. We begin by introducing the Lebesgue spaces 

(2.1) L2(17) ^ U- I \<l>{x)\^ da; < cx)| and Loo(^) = U ■ sup \(j){x)\ < ooj , 

and the Sobolev and Hilbert spaces 

(2.2) W^(f7) = {</)£ Lp{n) : D°(/) € Lp{n), for |a| < k} and H'=(f^) :== W^(n). 
These are equipped with the norms 

(2.3) ml^n^ - / 101' da;, Uk^^n) = ™P W^)\ . 

(2-4) Ikfwjcn) - E l|D"-rL,(a) ^^d l^^^.^^^ = J^ IP°-IIL(o) • 

|a|<fc \a\=k 

where o: = {«!,..., a^} is a multi-index, |a| = ^j^j^ a^ and derivatives D° are 
understood in a weak sense. We pay particular attention to the cases A: = 1, 2 and 

(2.5) Hj(f7) := closure of C;f (f^) in Hi(f]). 

The model problem in strong form is: Find u e H^(J7) n Hj(rj) such that 

(2.6) (^7., 0) = (/,(/)) V0eHj(17), 

where the data / G L2(r2) is prescribed and ^ is a general linear, second order, 
uniformly elliptic partial differential operator. Let A g Loo(f^)'^^'^, we then define 

^ ' u ^ ^u -.^ -A: B^u. 

We assume that A is uniformly positive definite, i.e., there exists a 7 > such 
that for all x 

(2.8) y'A{x)y>j\y\^ VyeR'*, 

and we call 7 the ellipticity constant. 
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Nonvariational PDEs are not as well studied as their variational brethren from 
a numerical analysis viewpoint. For the benefit of the reader we present a concise 
account of known results for strong solutions of this class of problem. 



2.1. Definition (strong solution). A strong solution of (1.1) is a function u £ 
H^(r2) n Hj(r2), that is a twice weakly differentiable function, which satisfies the 
problem almost everywhere. 



2.2. Theorem (existence and regularity of a strong solution of (1.1) |GT83} Thm 
9.15]). Let n C R"^ be a &'^ domain. Suppose also that A e C^inY""^ and 
f G L2{fl) such that the problem 

-A: J:>^u = f in Q 
(2.9) 

u = on ofl 



is uniformly elliptic. Then (2.9) has a unique strong solution. There also exists a 
constant independent ofu such that 

(2.10) ll«llH^(o)<C||/|lL,(a)- 

2.3. Remark (less regular solutions). Note that the theory of viscosity solutions 



has been developed for non classical solutions of (2.6) if the problem data does not 
satisfy the regularity assumed above see |GT83| . 

2.4. Assumption (regularity of A). From hereon in we will assume that the prob- 
lem data is sufficiently smooth such that solutions exist and belong to at least 

H2(f7)nHj(f7). 



2.5. Remark (regularity of fi). Theorem 2.2 specifies that Q must be a C ' do- 



main. We will be approximating such a domain with one which is only C ' (i.e., a 
polyhedral one). We thus assume that the model problem admits a unique strong 
solution even when f2 is only C ' . To circumvent this assumption curved finite 
elements could be used to fit the boundary exactly |Ber89| . For simplicity we will 
not present this case here, although we believe our analysis can be extended to this 
case. 

3. Discretisation 

Let =^ be a conforming, shape regular triangulation of O, namely, cf?' is a finite 
family of sets such that 

(1) K G ^ implies K is an open simplex (segment for rf = 1, triangle for d = 2, 
tetrahedron for d ~ 3), 

(2) for any K,J £ £^ we have that K O J is a full subsimplex (i.e., it is either 
0, a vertex, an edge, a face, or the whole of K and J) of both K and J and 

(3) [jKe^K = n. 

We use the convention where h : Q ^ R denotes the meshsize function of ,^, i.e., 
(3.1) h{x) :— maxhK, 

K3x 

where Hk is the diameter of K. We let S" be the skeleton (set of common interfaces) 
of the triangulation ^ and say e € (f if e is on the interior of fl and e G dfl if e lies 
on the boundary d^. The assumptions on the tessellation made here are typical 
in the finite element analysis. For the presentation of the method and its analysis, 
some assumption could be relaxed (e.g. the form of the elements or the assumption 
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on a conforming triangulation) but this would lead to an unnecessary increase in 
the complexity of the presentation. 

Let P'^{^) denote the space of piecewise polynomials of degree k over the tri- 
angulation c^ji.e., 

(3.2) P\^)^{(j,: ^UePHK)} 
and introduce the finite element spaces 

(3.3) Vz5 = Vi3(^,fc):=P'=(^)nHj(r!) 

(3.4) Vi5=Vz5(^,fc):=P^'(>^) 

to be the usual spaces of discontinuous piecewise polynomial functions which are 
compactly and non compactly supported respectively. 

3.1. Remark (generalised Hessian). Assume a function v e H^(il), let n : dfl — > 
R"^ be the outward pointing normal of fl then the Hessian D^u of v, satisfies the 
following identity: 

(3.5) I H^v (j)Ax = - I Wv®W(j)dx+ I Wv®n<l)ds V^eH^f]). 
JQ. Jn Jdn 



If V e H (Q) (3.5 1 is still well defined in view of duality, in this case we set 



(3.6) (D^i;!^) = - / Vu®V0da;+ / Vi;®n</>ds V G H^(rj), 

where the last term is understood as a pairing between H~^"(ri) and H^"(51). 

3.2. Definition (broken Sobolev spaces, trace spaces). We introduce the broken 
Sobolev space 

(3.7) }l''{£r) := |(/, : 0|^ g }i''{K), for each K e £^} . 

We also make use of functions defined in these broken spaces restricted to the 
skeleton of the triangulation. This requires an appropriate trace space 

(3.8) r(^):= n L2(9i^)= n H'(^). 

3.3. Definition (jumps, averages and tensor jumps). We define average, jump and 
tensor jump operators for arbitrary scalar functions v € T{S'), vectors v € T{S') 
and matrices V £ T(^) as 

(3.9) iv}^liv\K^+v\K,), iv}=^{v\K,+v\K,), 

(3.10) Ivj ^v\K^nK^+v\K^nK^, Ivj ^{v\k J'' uk^ +{v\K^)''nK^, 

(3.11) lV} = V\KinK,+V\K2nK2, H^ = 'wIki » "Xi + tll^a ® "^Ta- 
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Note that on the boundary of the domain dQ the jump and average operators 
are defined as 



(3.12) 
(3.13) 
(3.14) 



vj 



dn 



Iv} 



on 



■.= v, 



H 



m 



\an 



dn 



Vn, 



dn 



an 



We will often use the following Proposition which we state in full for clarity but 
whose proof is merely using the identities in Definition |3.3[ 



3.4. Proposition (elementwise integration). For a generic vector valued function 
p and scalar valued function (p we have 

(3.15) Y] / div{p)<j)dx= V(-/ p'^Vh(t>dx+ I 4)p'^nKds], 

k^^Jk ^,g^V Jk J ax ^ 

where V/i = (D/i)^ is the elementwise spatial gradient. Furthermore, If we have 

p e T{S U dVlf and 4> G T{S U dVl), the following identity holds 

(3.16) 

V / cj^p^UK ds = / bl I 0} ds + / W Ip} ds= f IpcPl ds, 
j^^g-Jaii Jg Jsuan Jsvjan 

An equivalent tensor formulation of ^3.15 )-(3.16) is 

(3.17) y^ / ^hP4'dx= V" (- / p<8)Vh0da?+ / (/)p (g) n/^ d^ 

r.-^ or Jk 1^^ OT^ JK JdK 



where 

(3.18) 

y^ / (jyp^riK ds 
KdSrJSK 



K^S 



Me 



ds- 



p1 ds 



gLtaii. 



iP4 



ds. 



suan 



In addition for matrix valued V we have that 

(3.19) V / [Bhp): Vdx= V (- / p: BhV dx + f {Vpynds 

Ke^rJ^ R-^ar^ JK Jan 



Kesr 



and 
(3.20) 



^ / [VpYnds = / IVf ip} ds+ / bl^: I V}ds= f IVpj ds. 
i^^^Jan Jg Jguan Jguan 

3.5. Construction of an appropriate discrete Hessian. We now use the frame- 
work set out in [ABCM02 to construct a general notion of discrete Hessian. We 
first give a definition using a flux formulation: 

3.6. Definition (generalised finite element Hessian: flux formulation). Let u e 
H2(5^) and U : H\5^) -^ Tits' U dfl) be a linear form and p : ^(5^) x H^(^)'' -^ 
T{(S U dQ) a bilinear form representing approximations to u and Vu over the 
skeleton of the triangulation. Then we define the generalized finite element Hessian 
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H[u] as the solution of 

(3.21) I H[u]<^dx = - j p®\/h<^Ax+ j p^(8)n$ds V$eHi(5^) 

JK Jk JdK 

(3.22) p<»qdx = - uBhqdx+ q (^ n Uk ds V q G (H^(^))'* , 
Jk Jk JdK 

for all $ e V_D. 

We now present the primal formulation for the generalized finite element Hessian: 

3.7. Theorem (generalised finite element Hess ian: primal form). Let u e y?{S^) 
and let U and p be defined as in DeGnition 3.6 then tie generalised Enite element 
Hessian H[u] is given for each $ e ^d as 

(3.23) 
[ H[u] $da; = - [ \/hU^Vh'^dx+ [ m^lPJ ds + [ i<S>im<g, ds 

Jfl Jq. J gudil Js 



«L/-u}[V/.$Lds 



svjdn 



U-u 



^h^J ds. 



Proof Note that in view of Definition |3 .3| for generic vector fields q G W and u G V 
we have the following identity 



(3.24) V / vq®nds= I H«)|g}ds 

K£sr-'9K Jsudn 



vJH^ds. 



Then summing (3.21 ) over K ^ £^ and making use of the identity (3.24) we see 
(3.25) 

/ H[u] $ da; = V" / H[u] $da;= V'f-/ P® V,i<i> dx + p^ (g) n $ ] di 
Jn j^^ or Jk r-c ^ ^ Jk J ok ^ 

m^lVK^ds^ \ m\VK\<^ds. 
Jg 



Kesr 



Kesr 



= pi^Vh^dx - 

'n Jgudn 



Using the same argument for (3.22) 
(3.26) 



/ p® qdx — 2_, / P'S) qdx — ^J f — / u D^q dx + q®n Uk dt 

Jn ucz-^Jk ircz^^ Jk JdK 



Kes 



Kesr 



= - / uBhqdx+ / [/ ®«g:& d,s+ / iUlM^ds. 

Jn Jgyjdil "- -" Jg 

Note that, again making use of (3.24) we have for each q G H^(5^)'^ and v G 

Hi(^) that 

(3.27) f q®Vhvdx = - j T>hqvdx+ j I q^ ^M ds+ f M^iv^ ds. 
Jn Jn Jgudn Jg 

Taking t; = w in (3.27) and substituting into (3.22) we see 

(3.28) 



pi^qdx= / qi^Vhudx 
n Jn Jgudn 



U-u 



iq}ds+ {{U-uJM^ds. 
Jg 



Now choosing q = V^.^ and substituting (3.28) into (3.21) we arrive at the fully 
generahsed finite element Hessian given by (3.23). D 
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3.8. Remark (consistent representations of the gradient operator). If one were 
interested in consistent representations of other derivatives, for example the gra- 
dient operator, one would need to modify the proof of Theorem |3.7| Exam- 
ples of consistent gradient representations can be found in |ABCM02] . See also 
|BO09[ iDPElOl IBE08| . Using this methodology it should be possible to construct 
an entire hierarchy of derivatives. 

3.9. Example. An example of a dG formulation for the approximation to the Hes- 
sian, D^u, can be derived by taking the fluxes in the following way 

(3.29) ^ rn ".lover ^ 

1^0 on 9il 

(3.30) P=iyhUh} onr^uan, 

where 9 G {—1, !}• The result is a discrete representation of the Hessian H[uh] as 
unique element of V^^ "^ such that 

/ H[uh] <Pdx = - / VhUh'^Vk'^dx 
Jn Jn 

)iVh'^}+m®i^hUh} ds 

Jn Jg 

)| V,,$} ds V$G Vc. 



(3.31) 



isudn 
3.10. The discontinuous nonvariational finite element method. We are now 



to find Uh e ^D together with H[uh] e V^^"* such that 



in a position to state the numerical method for the approximation of ( 1.1 ). We look 



(3.32) s^hiuh.-^) ^ l{-^) V*eVc 
with 

(3.33) ^h{uh,^) := f -A: H[uh]^'dx + f ah-'luhf m ds 

Jn Jsuon 

(3.34) l{^):= f f^dx, 

Jn 

where the penalisation parameter cr > is to be chosen sufficiently large to guar- 
antee coercivity. 

Using the L2 projection operator Py : L2(il) — > ^d defined for v e L2(fi) through 

(3.35) / Pv(w) ^ da; = / w* da; V* G V_d 
Jn Jn 

it is possible to elliminate the finite element Hessian from the bilinear form for 
sufficiently smooth A: 
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3.11. Lemma (elimination of the finite element Hessian in a g eneral setting). If 



1 dxd 



r k■4-^ 1 
A e W;^ (ri) and the Suxes are chosen as m Example 

(3.36) 



3.9 



then 



^h{uh. *) - / Dh{Vy{-^A)) VhUh dx- f ^Kf I ^/.(Pv(*A)) } ds 

Jsuaa 



'Sudn Jguon 

Proof. This follows from the following identity 
(3.37) 

f -A: H[uh]'i'dx= [ -H[uh]:{'i'A)dx= [ -H[uh]:Pvi'i'A)dx 
Jn Jn Jn 

= / DhiFvi^A)) VhUh dx- j eiuhf I DH{Vy{-^A)) } ds 

-f lPv{^A)f iWhUh} ds. 
J Sudn 

D 

3.12. Remark. The solution of the problem in this form is nontrivial due to the 
global L2(ri) projection appearing in the formulation. However, in the discontinuous 
setting the global L2(il) projection is in fact computable locally. We may actually 
exploit this fact to optimise our schemes efficiency. We will discuss this further in 
the sequel. 



3.13. Example (Laplacian formulation). Note that if in (1.1) we have that A ■ 
then we have that 

(3.38) / = -^:D2u = -Au 

and our bilinear form reduces to 

s^h{uh,^)= j {Vh-^yVhUudx- f 0Kf {{V,,*} ds 
/o oq\ Jn J Sudn 



mi^hUhV -<yh-^iuhVmd 



s 



isuon 
since Pv(*A) = */. 

The nonvariational finite element method thus coincides with the classical (sym- 
metric) interior penalty method for the Laplacian |DD76j . 

3.14. Remark (relation to standard dG methods). It is not difficult to prove that 
choosing to numerical fluxes in the same way as presented in |ABCM02l Table 3.2] 
results in the same correlation to the dG methods summarised in the aforementioned 
paper for the case that A is constant. For brevity we will not prove this here. 

Note that when A is not constant we have that the nonvariational finite element 
method does not coincide with its standard variational finite element counterpart. 
There is an extra stability property which allows the method to successfully cope 
with advection dominated problems |LP11[ §4.2] which is illustrated by the result 
of Lemma 13.111 

We conclude this section with a proof consistency of the method and then show 
that Galerkin orthogonality holds. 
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3.15. Lemma (consistency). Let m G H {,'7) and assume that the numerical Buxes 
are chosen in a consistent fashion in the sense of |ABCM021 §3.1], that is, 

(3.40) U = u\sudn 

(3.41) p = Vwl^uao 
Then for $ G V^, 

(3.42) / H[u]^Ax^ / D^-u$da; 

Therefore we have that H[u] = Pv(D^u). 

Proof Applying Proposition |3.4| to the first term in the definition of H[u] yields 

(3.43) 

/ H[u]<^dx= I D2ii$da;+ / Ip-Vu]^|$} ds + j |p- Vu}«)|$l ds 

Jii Jn Jg Jsudn 



f lu~u}l\/^^^ds- I 
Jg Jg 



B^u^dx V$eV 



gudn 

D- 



U -u 



V,,$| ds 



which proves the results under the consistency conditions on the fluxes. D 

3.16. Lemma (Galerki n or thogonality). Let u G H'^(il) n Hj(r2) he a strong so- 

r I o 

lution to the problem {1.1) and let Uh G V^i he its nonvariational finite element 

approximation. Assume that the numerical fluxes U and p are consistent then we 

have the following orthogonality result: 

(3.44) £/h{uh-u,^) = J{^) V*gVz5, 
with the error functional given by 

(3.45) J(*) = / {D^u - H[u\) : (A*) da;. 

Jn 

Proof Using the consistency result and that Jm] = we conclude 

£/h{uh - u, *) = ^^(u,,, ^) + I A: Hlul'^dx = Z(*) + / H[u] : (A^) dx 

Jn Jn 

= - A: D^u^ - H[u]: (A'i) dx = J(^), 
Jn 

concluding the proof. D 

3.17. Remark. If A is piecewise constant then since H[u] = Pv(D^u) we have 
J(^) = and we recover the usual Galerkin orthogonality £4(u/i — u, ^) = 0. We 
will show in the next section that in general the error functional J is of higher order 
for smooth enough u. 

4. COERCIVITY, CONTINUITY AND CONVERGENCE 

In this section we examine the coercivity, continuity and convergence of the 
method. We will focus on the fluxes given in Example |3.9| to simplify the presen- 
tation. Furthermore we make the following additional assumption on the problem 
data. 
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4.1. Assumption (coercivity of the nonvariational problem). For the rest of this 
section we will assume that the nonvariational operator is coercive, that is A G 

[W^^(r2)]'""^ and that div(DA) < 0. 



4.2. Remark (variational nature of the coercive problem). Under Assumption 4.1 



the problem can be written variationally. The solution to the nonvariational prob- 
lem is the minimiser to the (degenerate) second order variational problem: Find 
u e H^(n) n Ho(f]) such that 

(4.1) ^M= inf /H, 

where 

(4.2) ^[v]:^ [ {A-.D'^u + f-D^: A)udx. 

4.3. Definition (H^(3^), H^(5^) and H"^(^) norms). We introduce the broken 
H^^) and H^(,5^) norms as 

(4-3) IWhfdG,! — W^ hUhWt^^n) + f^^^ Wl^hlWl^^s) ^ 



and the H ^{3^) norm as 

(4-5) \\uh\\dG-i ■= sup 



Jj-j UhVh dx 



These are equivalent to their continuous equivalent norms for functions in V/j. 

4.4. Proposition (projection approximation in Vd). Let Py : L2(^) -^ ^d be the 
L2(rj) orthogonal projection operator defined by (3.35). Using standard approxi- 
mation arguments we have that 

||w - Pv wllrfG^i <C/i'= |w|Hfc+i(!:2) and 

ll« - Pv «IIl2(o) < C^fT'''^^ blH'=+i(a) ■ 

In particular, let Ah denote the L2 orthogonal projection of A into the space of 
piecewise constant functions, then we have 

(4.7) II A - Ah\\^^^^^^ < (1 + C^h) ||DA||L^(f,) . 

4.5. Theorem (stability of H |Pryl2[ Theorem 4.10]). Let H be defined as in 
Example \3.9\ then the dG Hessian is stable in the sense that 

(4.8) plvh - ffK]||L,(o) ^ ^(/ ^'^ llVhMl' + h-^\lvh\f ds 
Consequently we have 

(4.9) ll-ffNllL(o)<C|k'.llL.2- 
We now state the following technical Lemmata. 
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4.6. Lemma (upper bound on \\vhA - Py{vhA)\\^^f^^^^). Let A e [w^"^(ri)l 
and Vh &y then it holds that 

(4.10) 

WB^iv^A - Pv(z;,A))||l,(^) < C2hUl + C,h) ||DA||l^(^) + f] |A|^.^(^) j W^hvhW^^^^^ 

Proof Let Ah denote the L2 orthogonal projection of A into the space of piecewise 
constant functions. Then adding and subtracting appropriate terms we see 

(4.11) W^hivhA - Pv(«„A))||l^(j,) = \\Bn{vh{A ~ A^) - Py[vnA - A^M^^^^^ . 

Using the approximation properties of the L2(ri) projection given in Proposition 
14.41 we see that 



\Bh{vhA - Pv(i;hA))||L,(^) < Ch'' \vh{A - A,0|h.+i 



(O) 



(4.12) ^ 

<Ch II A - AhWy^k+i^^-^ hhWnk+^ii) 



Now using the properties of A^ (4.7) and inverse inequahties we have 

(4.13) 

\\BhivhA - Pv(fhA))||L^(j^) <Ch\\A- Aft||wfc+i(jj) W^^hVhh^^n) 

< ChUl + C,h) ||DA||l,^(,,) + ^ |A|w^(^ J ||V,t;,||L,(o) 

as required. D 

4.7. Lemma (upper bound on some skeletal terms). Let Vh G "^ d, B e W^(J7)''^'' 
and Bh & V^^ be the L2 orthogonal projection of B, then in view of trace and 
inverse inequalities we have the following: 

III DhPy{vhB)MLAg) < Ch-^'^ ||D„(Pv(z;,.S))||l,(^) 

<C3/i-i/2||B||w^(f,)||t.„LG,i. 
||IPv(«,B)1IIl,(<,) < Ch''^ ||D;,(Pv(«^B))||l,(^) 

(4-14) <C^h^/^\\B\\^^^^^^\\vhha,i^ 

III V/.«ailL,(^) < C^h-^'^ W^hVhWi^^i^n) > 

||« «ailL.(^) < Ge/i'/' llVhi'/.llL.co) : 
||ISJ||l^(^)<C7/i||DB||l^(,,). 

Proof For brevity we prove only the first inequality, the second and third follow 
similar arguments. In view of the definition of the average operator (3.91 it follows 
that 

(4.15) III D^ Pv(z;;.B) }||l,(^) ^\T. H^/^ Pv(«hB)||L,(9^) . 

Now by a trace inequality we see that 

(4.16) III D, Vy(vuB) :&||l,(^) < C ^ h-^'^ ||D, Pv(«/.B)||l,(k) ■ 

K£sr 
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Using the stability of the L2(ri) projection operator in H (fl) we have 

(4.17) Kes 

<C3/i-i/2||B||w^(o)ll«'.LG,i' 
as required. D 

4.8. Theorem (discrete continuity and coercivity). Let the conditions in Assump- 
tion \4.1\ hold. Suppose also that DA is sufficiently small such that 

(4.18) 
7 - 2e - (C4C5 + Cidh + ^^h + C2h{l + Cih) ) \\V)A\\^^^^^ 

-C2/iEl^lw.^(O)>0, 

4 = 2 



where 7 is the ellipticity constant, Ci is a constant appearing in Proposition 4.4 



with fc = 0, C2, C4 and C5 are the constants appearing in Lemmata \4. 6\ and \4. 7\ and 
e > is some parameter. In addition assume a the penalisation term is sufficiently 
large, specifically 

(4.19) 

C|(l + C,hf ||DA|1^^(^) + 4e2 / cc^il + C^h) 

'^ 4e 

(e + ifclwiA, 



[ ' "\ ^ ' + G,h)^ IID^IIl^(o) 



2 



>0. 



4e 

where Ah is some piecewise constant approximation to A. 
Then there exist positive constants Cb a-nd Cc such that 

(4-20) \-s^h{uh,Vh)\<CB\\uh\\^iQi\\vh\\aQ-^ and 

(4.21) ^h{uh,Uh)>Cc\\uh\\la,i Vu/,.,u;i e Vz,. 

We postpone the proof of this theorem to the end of this section and first prove 
the error estimates for our discrete solution. 

The results of Theorem 4.8 allow us to invoke Strang's second Lemma. 

4.9. Corollary (Strang |EG041 c.f.]). Tiiere exists a C > such that 

tAoo\ II II ^ n( ■ f II II , Wh{u,Wh)-l{wh)\\ 

(4.22) \\u-uh\\^a.i^'^\ inf \\u-vh\\^a,i+ ^up r—- . 

\vKeVo ni^eVn ll"'''ll<iG,l / 

4.10. Lemma (discrete negative norm convergence). Let A G Loo(i^) and u £ 
H "'''^(rj). Then we have that there exists a constant C > such that 

(4.23) II A: (D^u - /f M)Lg,-i < G^'+' Halloo IMh-^H^) ■ 



Proof We have, in view of Cauchy-Schwarz inequality, that for ^ e V^ 
(4.24) 



f A:{-D^u~H[u])^dx<\\A\\^^^,,^\\B^u-H[u]\\^^^^^mk^^^) 



<Ch'+'\\AL^in)\Mu''^H^)mUi^ 
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since H[u] = PyD^u and by the definition of H'II^jq i- The result follows noting the 
definition of the discrete negative norm in (4.5 1. D 

4.11. Remark. Noting the definition of the error functional in the Galerkin orthog- 
onality we deduce that 

(4.25) ||j||^^_^^0(/.^-+i). 

4.12. Theorem (convergence of the nonvariational method). Let u solve the non- 
variational problem {1.1) and {uh,H[uh]) solve the nonvariational Bnite element 
approximation {3.32) where H[uh] is a consistent approximation oflfiu (for exam- 
ple that given in Example 3.9). Then the following error bound holds: 

(4.26) \\u - u,,||^g_i < C{h'' |m|h'=+ho) + h''+' Mn^+^n)) ■ 

Proof The proof of (4.26 1 is immediate from applying Proposition 4.4 to Corollary 



4.9 with Vh = P\/u and noting that the bound for the consistency error nothing but 



the result of Lemma 4.10 concluding the proof. D 

To conclude this section we prove Theorem |4.8| 
Proof Theorem |4.8|Let u^^Vh E Vjj, then we have 



■^h{uh,Vh) ^ - A: H[uh]vhdx + ah ^ / lu/J"^ Ivhj ds 

Jn Jsuan 

(4.27) = - / H[uh] : {v^A) dx + ah'^ f lu^V KJ ds 

= - f H[uh] : PvivhA) dx + ah-' f lu^f Kl ds. 
Jn Jsudn 

Now making use of the definition of H from Example |3.9| we see 

(4.28) 

£^h{uh,Vh) = - I ^luh:Vy{vhA)dx+ / iV/^uJ^g : | Pv(v/iA) } ds 
Jn J g 

eiuhf I D;,(Pv(«„A)) } ds + cjh-' f Kf Kl ds. 

Jsudn 



I (Sudn Jsudn 

Adding and subtracting appropriate terms we have that 
(4.29) 

£/hiuh,Vh) ^ ~ / A:T)luhVh + T:)luh:{vhA-Py{vhA))dx 
Jn 

+ f iVhU J^ : I PvivhA) Ids^ f 9luhV | D^(Pv(i'„A)) } ds 
Jg Jgudn 



g Jgudn 

-1 



+ ah-' / luhf Kl ds, 

Jgudn 

which rewriting variationally gives 
(4.30) 



■i^hiuh,Vh) = - / BhiA'^hUfi) Vh + BhA'^hUhVh + D^u^ : {vhA - Pyl^VhA)) dx 
Jn 

+ f iVhW J^ : I PvivhA) Ids- f eiuhf I D^(Pv(«„A)) } ds 
Jg Jgudn 



g Jgudn 

-1 



'^h-' / lu,,r Kl ds 

Jgudn 
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Note that 
(4.31) 



- / BhiAVhUh) Vh + Bluh : KA - Pv{vhA)) dx 
Jn 

Y^ / {AVhUh)yhVh~'Dh{vhA-P^{vhA))VhUhdx 



Kesr 



+ / -{A^hUh)vhn+{{vhA-Vy{vhA))^hUh)nds 



I OK 



^ / {AVhUh)VhVh-^h{vhA-Vy{vhA))VhUhdx 



Ki^Sr 



IdK 



{¥y{vhA)SJhUh)nAs 



/ [AVhUh) ^hVh - Dh{vhA - Vy{vhA)) VhUh Ax 
Jn 

iPy {vhA)f iVhUh} ds~ f lVhUhl^:iPvivhA)}ds, 
I Jg 



gudn 



and hence we see that upon substituting (4.31) into (4.30) that 
(4.32) 

J^h{uh,vh) = / {AVhUh) VhVh + DAVhUkVh - BhivhA - Py{vhA))VhUh dx 
Jn 

eiuhf I T)h{Pv{vhA))} +lPy{vhA)f I WhUh} ds 



isvjdn 



gudn 



Kf Kl ds. 



We proceed by applying Cauchy-Schwartz componentwise to (4.32) and estimating 
by 1 arriving at 

(4.33) 

s^h{uh,Vh) < l|V,,u/i||L,(o)(M|JL^(n) l|V;,w/,||L3(t2) + ||DA||l^(o) II^/iIIl^ch)) 

+ WBhivhA - Py{vhA))\\^^^^-^ I|V/.u/.||l,(o) 

+ IIKlllL.(<f)ll{{D/.(Pv(z;/.A))}||L,(^) 

+ WlPWhAJW^^^g^ III VhUh}\\^^^g^+ah-^ IIK1IIl,(^) II WIIl.(^) • 
In view of Lemma |4.6| and the Poincare inequality we have 
(4.34) 
^h(uh,Vh) < (II A||l^(^) + Cp ||DA||l^(^) + C2h II A||^^+i(f,)) ||V,,w„||l^(j,) II Vhf/.|lL,(o) 



iud 



\L2{S) 



Dh{Py{vhAj) 



IIl2(.?) 



IIIPvWhAlllL,(<?) Ill "^ hUhJh^t^s) + <jh 1 ||K1IIl.(<?) IIKlllL2(,f) 

we apply the result of Lemma |4.7| wl 
e that 

\^hiuh,Vh)\ < Cb \\uh\\dG,l II^'^Lg,! 



For the skeletal terms we apply the result of Lemma |4.7| which upon substituting 
these into (4.34) we see that 

(4.35) 

as required. 
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For coercivity we use the equality given in (4.32 1 with Vh — Uh to find 

(4.36) 

J^h{uh,Uh) = I {AVhUh) ^hUh + T)AVhUhUh - D,j(u;,.A - Vy{uhA))VhUh da; 
Jn 

eiuhf I Bh{Vy{uhA))^ +l¥y{uhA)f I V,ua 



isuon 
+ ah-^luhV luhl As 

6 

i=l 

We proceed by bounding each term individually. By the ellipticity of the problem 
we have that 



(4.37) J^i= I {AVhUh)VhUhdx>-i\\VhUh\\l^^n)- 
By the coercivity of the problem we have 

(4.38) ^2=1 ^AVhUhUh - / V)A\Vh{uj) Ax ^ -\ [ div(DA)M^da; > 0. 

Jn Jn ^ ^ Jn 

By the Cauchy-Schwartz inequality and making use of Lemma |4.6| 



-Jz= \ V)h{uhA-Vy{uhA))^ hUhA'x. 
Jn 

(4,39) < \\Qh{uhA - Pv(u,,A))||L^(f^) l|Vh.Wh||L2(n) 

< C^hU + C,h) ||DA||l^(,,) + ^ |A|^^(,,J |lV.t-,|l^^(^) . 

We combine the fourth and fifth terms and let Ah denote the L2 orthogonal pro- 
jection of A onto the space of piecewise constant functions. Upon adding and 
subtracting appropriate terms 

(4.40) 

-^4 - ^5 - / ^Kf {{ ^h{^y{v^hA)) S +lPv(ii/, A)f I ^hUh } As 
Jsudn 

0{luhf I Dh{Pv{uhA - UhAh)) } +0Kf I ^hiuhAh) }) 
+ iPviuhA - UhAhW {{ "^hUhl +KAJT I VhUhJ As. 
Using the identities 

(4.41) / KA^l As^ f lAh} luhl As+ f IA,f fuh} As and 

Jsudn Js Jsudn 

(4.42) 

/ I AhVhUh } As^ f ] lAhl iVhUhl As+ f | A^ }{{ VhUh } As 
Jsudn Js 4 Jsuon 
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we have that 
(4.43) 

Jsudn 
+ lVy{uhA-UhAh)f I Vhunl +(^+ 1)(| A,} KJ) | V,,..,,} 

Using Caiichy-Schwartz we see 

+ \\lVy{uhA - UhAhM^^i^g-^WlV hUh1\\u(g) 

(4.44) +(^ + 1) III ^aiiL^(^) iikiiIl.(^) III ^hUhiWuis) 

+ iii^ftiiiL„(^) III ^hUhih^ig) III ^iaiiL.(^) 

+ 4 IIIAlllw,(<?) IIKIIIl.(^) IIIVhwJ|lL,(^) . 

Making use of the various bounds from Lemma |4.7| we have 

(4.45) 

-^4 - ^5 < C73/1-1/2 p _ Aft.||^i^(^) h/^Lci IIK1IIl.(,?) 

+{o+i)c,h-^'^ III AaiiL^(^) iikiiIl.k) iiVh^^^iiL,(.i) 



C'sC'eC'y/l 11 11 II ||2 



D^llL^m) l|Vh"/JIL(m +^7 l|DA||, .^. ||K1| 



2 



We now apply a Cauchy inequahty and use the approximation properties of Ah 
from Proposition 4.4 to find for any e > that 

JTi = Cs/i"^/^ ||A - Ah\\^fgi^(^) \Wh\\dG.i IIK.1IIl2(,?) 

C3 ll^^^/illwj^cn) ,_i 



< 



4e 
(4-46) ^C?||^-^,|k^„wo.+4e 



h IIK1IIl2(<?) + *=II"'»Lg,i 



^r^ h IIKlllL,(^-) + e||v,jWft|| 



< 



C|(l + Ci/i)'||DA||^ m)+4e' , , , 

^^^ /^-MlKlllL(^)+^l|V/."/JlL(a) 



4e 
The other terms are bounded similarly in that 

JC2 = CiCr^ \\A - ^h|lwj^(s:2) ll^'JIdci ll^'i'"'illL2(n) 
QCsP-A/ill^i^^j^) 

< ^ (ll^'JIdci + W^hUhU^^n), 

(4.47) ^ , 

< ^4^5(1 + Ci/J) ||DA||l^(^) II V,K;,||^^(^) 

+ ^ h ||K1IIl2(^)^ 
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and 



(4.48) 



=^3 =(e + 1)^5/^-'/' ii{{ AaiiL^(^) iikiiIl,(^) iiVhWhi 



< e\\WhUh 






lL2(^) 



lL2(n) 



'"iiL2(!:2) -r ^ 

Note that the final two terms are already in their desired form since 



/i"MlKll!L(^) 



(4.49) 
(4.50) 



Jza — 



l^^k^m W^hUnWl 



2(a) 



JTs = C7 \\DA\ 



Loo(fi) 



2 

L2(<f) 



Collecting the bounds from (4.46)-(4.50) shows 

(4.51) 

CeCi 



^4 - A < 2e + Cs C4 + /i C1C4 + 



IDA! 



L=o(0) 



V/iU/il 



L2(0) 



'C|(l + Ci/i)^||DA||^^(,,)+46 



4e 

6-4^5(1 + C1/1) 



C7/1 IIDAI 



Loo(f2) 



1)"C| 111^,411^ m\ , 

^^^^U-MIK1I1l2(^) 



4e 



The final term in (4.36) is given by 

(4.52) 



^6 = / ah-'luhf Kl ds = a/i-i |1K1|1 



Finally, collecting the bounds from (4.37), (4.38), (4.39), (4.51) and (4.52) shows 



(4.53) 



fe+i 



^hiuh,u>,) > 7 - 26 - C2/i (1 + C^h) ||DA||l^(,,) + ^ \A\ 



w^(n) 



C5{Ci + h)(ciCi + 



^ ^ "DA||L^(f2) ) IIV/.u/JIl^^jj) 



Cl{l + C,hf\\BA\\l,,,.+Ae 



4e 
6*4^5(1 + Ci/i) 



Ct/i ||DA|i 



Loo(n) 

4, ,'^"'IIK1IIl2 



i)ViiiSAaiiL (^)\ 1 2 



Coercivity of the discrete bilinear form follows using the assumption in Theorem 
4.8[ by choosing e sufficiently small and the penalisation parameter a sufficiently 
large for small enough h. D 
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4.13. Remark (the coercivity bound). We note that the coercivity bound rehes on 
the term DA not becoming too large, as specified in Theorem |4.8[ If it is we view 
this as an advection dominated problem. Our numerical experiments suggest that 
there is no condition on the size of this term. 

If the coefficient matrix A is divergence free, i.e., DA = then the bound 
simplifies considerably. For example, in the case that A is constant we regain the 
same theoretical results as for the method given in Example |3.13[ 

5. Numerical experiments 

In this section we detail numerical experiments carried out in the finite element 
package Du ne-Fem |DKN010| which is based on the Dune software framework 
|BBD+08a[ lBBD+08b| . The code used to test the method will be made freely 
available within the DuNE-FEM-Howto in a future release. 

We present some benchmark problems designed such that the exact solution is 
known. In each of the experiments the domain ft = [0, 1]^ and we consider the 
coefficient matrix to be 

1 b{xy 
b{x) a{x)_ 



(5.1) Aix) 



varying a{x) and b{x). 

In each of the numerical experiments we make use a stabilised conjugate gradient 
solver taken from the Dune-Istl module }BB07j preconditioned with an incomplete 
LU factorisation. We choose the penalty parameter a = 20. 

5.1. Test 1 : a coercive operator. In this test we take the components of A 
such that the operator is coercive, fitting into the analytical framework presented 
in El With X —{xi,X2), we set 

(5.2) a{x) = - ln((a;i - 1/2)^ + 10"^") + 1 

(5.3) b{x)=0. 

We choose the problem data such that the exact solution is given by 

(5.4) u{x) ~ sin (nxi) sin (71x2) 



and approximate this using the formulation (3.32 1. In Tables la- lb we present the 
results for the cases k — 1,2, numerically demonstrating that the analytical rates 
of convergence are achieved in the dG energy norm, moreover, optimal convergence 
is achieved in L2(r2). 

5.2. Test 2 : nondifferentiable operator |LP11[ §4.4]. In this test we take A 
such that it is comparible to [LPlli §4.4]. We take 

(5.5) aix) = 2 

(5.6) b{x)^{xlxlf\ 



We choose the exact solution as in l5.1l and conduct the same tests. Tables [2aH2bl 
detail the results. Note that this is not a coercive operator and as such, does not fit 
into the analytical framework presented in fJ4J we do however still achieve optimal 
convergence in ||-|| and IMI^g i- 
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5.1 



Test 1. We present errors and convergence rates 



Table 1. 

of the approximation given by solving (3.321 



(a) Piece wise linears, A: = 1. 



^elements 


\\u— Uh\\ 


EOC 


IIw-"'>Lg,i 


EOC 


128 


0.0196123 


1.86116 


0.414643 


0.953491 


512 


0.00506166 


1.95408 


0.209225 


0.986817 


2048 


0.00128044 


1.98298 


0.104907 


0.995937 


8192 


0.000321803 


1.99238 


0.0525047 


0.998597 


32768 


8.06862e-05 


1.99578 


0.0262623 


0.999456 



(b) Piecewise quadratics, k = 2. 



^elements 
128 
512 
2048 
8192 
32768 



\\u — Uh\\ 

0.000475513 
5.99935e-05 
7.52887e-06 
9.42737e-07 
1.17929e-07 



EOC 
2.96865 
2.98661 
2.9943 
2.99751 
2.99893 



0.0308463 
0.00779373 
0.00195531 
0.000489487 
0.000122443 



EOC 
1.95 
1.98471 
1.99492 
1.99805 
1.99916 



5.2 - Test 2. We present errors and convergence rates 



Table 2. 

of the approximation given by solving (3.321 



(a) Piecewise linears, fc = 1. 



^elements 


\\u — Uh\\ 


EOC 


h-^hW^ci 


EOC 


128 


0.0172648 


1.89433 


0.41799 


0.955709 


512 


0.00441656 


1.96684 


0.210818 


0.987469 


2048 


0.00111269 


1.98887 


0.105688 


0.996186 


8192 


0.000278969 


1.99588 


0.0528915 


0.998707 


32768 


6.98234e-05 


1.99832 


0.0264548 


0.999507 



(b) Piecewise quadratics, k = 2. 



^elements 
l28 
512 
2048 
8192 
32768 



\\u — Uh\\ 

0.00047216 
5.98325e-05 
7.52118e-06 
9.42426e-07 
1.17933e-07 



EOC 
2.9534 
2.98028 
2.9919 
2.99651 
2.99841 



0.0309416 
0.00782197 
0.00196325 
0.000491575 
0.000122975 



EOC 

1.9514 

1.98394 

1.99429 

1.99776 

1.99904 



5.3. Test 3 : irregular solutions. In this test we consider the case the exact 
solution does not satisfy the regularity requirements presented in the analytical 
framework of MJ i.e., u ^ H """^(ri). In addition we consider the case that u ^ 11^(17), 
demonstrating the method converges even for viscocity solutions of the problem. 

We consider the coercive operator from §5.1| and choose the problem data such 
that 



(5.7) 



u{x) 



cos ( Stt a; 



1|2 
2 I 



1 



\X — 



if ^ 

otherwise . 



< 



DG METHODS FOR NONVARIATIONAL PROBLEMS 



21 



Note that this function is H^(f7) but not H'^(i7). We also take the problem data 

100xi(l-a;i)a;2(l-a;2) 



such that 
(5.8) 



u{x) = 



This function is H"'^(51) but not H^(f7). The results are given in Tables 



3a 



3b 



In the case u is given by (5.7) the scheme converges with optimal rate in the 
ll'lldG 1 norm even if the solution is not in H . The convergence in the L2 is more 
erratic, but we observe the same behavior testing the standard IP FEM taking A 
to be the identity. 



In the case u is given by (5.8) the convergence rates are suboptimal since the 
solution is not H^. 



Table 3. 



5.3 



Test 3. We present errors and convergence rates 

In both cases we 



of the approximation given by solving (3.32) 
consider k = 1. 



(a) The solution here is given in |5.7| . The function 
ne H2(Q) but u ^H3(n). 



^elements 


\\u — Uh\\ 


EOC 


ll"-'"h|ldG.l 


EOC 


128 


0.0362651 


2.47943 


0.837082 


0.939619 


512 


0.0267684 


0.43805 


0.406003 


1.04388 


2048 


0.0179914 


0.573227 


0.253977 


0.67679 


8192 


0.00292357 


2.6215 


0.103168 


1.29971 


32768 


0.00174473 


0.744729 


0.0541648 


0.929566 


131072 


0.000421749 


2.04854 


0.0258935 


1.06476 



(b) The solution here is given in I 5.8 I. The function 
ue Hi(Q) but n0H2(n). 



^elements 


\\u — Uh\\ 


EOC 


ll^-^hlldG.l 


EOC 


128 


0.223469 


1.80378 


6.42181 


0.843123 


512 


0.0616572 


1.85773 


3.49469 


0.877816 


2048 


0.017159 


1.84531 


1.87984 


0.894556 


8192 


0.00509901 


1.75067 


1.00295 


0.906363 


32768 


0.00177874 


1.51936 


0.531521 


0.916047 


131072 


0.00076433 


1.21859 


0.280092 


0.924224 



6. Conclusions and outlook 



In this work we have extended the framework from |LPllj for linear nonvaria- 
tional problems to incorporate discontinuous approximations. 

We have shown the method presented (and subsequently that of the continuous 
case from |LP11) ) is well posed and converges optimally under coercivity assump- 
tions on the coefficient matrix A. 

In the numerical experiments we note the the method is well posed and converges 
optimally even for A which do not satisfy the coercivity assumptions or u which 
do not satisfy the regularity needed in the analytical framework. This motivates 
that another analytical approach needs to be developed. This approach can not be 
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variational in nature as such will be completely non standard. This is the topic of 
ongoing research. 
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